Trait analysis for no competition pots
For the overall trait analysis, we are going to do two separate analyses, one for each data set. Because there are so many potential interactions between factors (environmental or genetic), we will use a backward selection model approach to see what interactions are important in explaining variation in the data. Here is the general approach:
- Fit the most complex model with respect to fixed effects
- 5-way interaction between \(\text{CO}_2\) (2 levels: ambient, elevated), salinity (2 levels: GCREW, freshwater), elevation (continuous), age cohort (2 levels: ancestral, descendant), and provenance (2 levels: Kirkpatrick Marsh, Corn Island)
- quadratic term for elevation (no interactions)
- initial weight of the propagule, temporal planting group (2 levels: late June, early July), lab of origin (2 levels: ND, UTK)
- In this model also include random intercepts for genotype and for chamber
- For some models, the variance explained by chamber is basically zero and can be dropped from the model
- Use
lmerTest::step for perform backwards stepwise selection on the fixed effects of that model
- Use the argument
keep to force the stepwise selection procedure to keep the additive fixed effects of \(\text{CO}_2\), salinity, elevation, age cohort, and provenance because these are things that we purposefully manipulated
- This means that the only terms that can be dropped are interactions and covariates related to the experimental set-up such as initial weight of the propagule
- Given the final model, then test to see if there are genotype by environment interactions that can explain significantly more variation
- Do this by fitting models with random slopes using terms (up to two-way interactions) that were kept in the model for \(\text{CO}_2\), salinity, and elevation
- Test to see if the genotype by environment should be kept using a significance test with
anova() while setting the argument refit = FALSE to fit models using restricted maximum likelihood for these comparsions
Aboveground biomass
# Center and scale continuous variables and subset to just aboveground biomass
# greater than 0
bg_nocomp %>%
mutate(elevation_sc = scale(elevation)[,1],
weight_init_sc = scale(weight_init)[,1]) %>%
filter(extinct == 0)-> bg_nocomp
# Fit the most complex model -- sqrt of agb biomass seems like the best
# transformation
agb_mod <- lmer(sqrt(agb_scam) ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|site_frame) + (1|genotype), data = bg_nocomp)
## boundary (singular) fit: see ?isSingular
# Warning indicates that one of the random intercepts is estimated to be zero
VarCorr(agb_mod) # It's chamber, so we'll drop it
## Groups Name Std.Dev.
## genotype (Intercept) 0.24224
## site_frame (Intercept) 0.00000
## Residual 0.51650
# Fit a model without that term
agb_mod2 <- update(agb_mod, .~.-(1|site_frame))
# Do a stepwise model selection
# step(agb_mod2, keep = c("co2", "elevation_sc", "salinity", "age", "location"))
# There's a model fitting issue because the model is too complicated. Let's see
# if there's a significant 5-way interaction
Anova(agb_mod2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: sqrt(agb_scam)
## Chisq Df Pr(>Chisq)
## weight_init_sc 4.5169 1 0.033562 *
## date_planted_grp 1.1259 1 0.288661
## origin_lab 0.0399 1 0.841596
## age 0.4286 1 0.512701
## location 0.2874 1 0.591898
## co2 4.2801 1 0.038561 *
## salinity 0.3342 1 0.563194
## elevation_sc 68.4369 1 < 2.2e-16 ***
## I(elevation_sc^2) 25.5516 1 4.307e-07 ***
## age:location 1.5009 1 0.220539
## age:co2 0.0240 1 0.876871
## age:salinity 0.1001 1 0.751663
## age:elevation_sc 0.0282 1 0.866728
## location:co2 0.6727 1 0.412121
## location:salinity 0.0117 1 0.913832
## location:elevation_sc 0.5459 1 0.459995
## co2:salinity 0.7830 1 0.376220
## co2:elevation_sc 4.1267 1 0.042212 *
## salinity:elevation_sc 0.0245 1 0.875622
## age:location:co2 2.4417 1 0.118148
## age:location:salinity 0.0659 1 0.797354
## age:location:elevation_sc 0.0807 1 0.776346
## age:co2:salinity 4.9315 1 0.026371 *
## age:co2:elevation_sc 0.0013 1 0.971258
## age:salinity:elevation_sc 0.0832 1 0.772974
## location:co2:salinity 0.5042 1 0.477659
## location:co2:elevation_sc 0.1599 1 0.689287
## location:salinity:elevation_sc 1.3568 1 0.244085
## co2:salinity:elevation_sc 0.0737 1 0.785970
## age:location:co2:salinity 0.3398 1 0.559959
## age:location:co2:elevation_sc 0.2862 1 0.592665
## age:location:salinity:elevation_sc 0.4013 1 0.526431
## age:co2:salinity:elevation_sc 7.6946 1 0.005539 **
## location:co2:salinity:elevation_sc 2.2696 1 0.131935
## age:location:co2:salinity:elevation_sc 0.0104 1 0.918683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Nope! Drop it and try stepwise again
agb_mod3 <- lmer(sqrt(agb_scam) ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^4 + I(elevation_sc^2) +
(1|genotype), data = bg_nocomp)
step(agb_mod3, keep = c("co2", "elevation_sc", "salinity", "age", "location"))
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 37 -207.59 489.19
## (1 | genotype) 0 36 -212.58 497.17 9.9812 1 0.001582 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## origin_lab 1 0.0103 0.0103 1 31.802
## age:location:co2:elevation_sc 2 0.0738 0.0738 1 183.667
## age:location:co2:salinity 3 0.0958 0.0958 1 177.380
## age:location:salinity:elevation_sc 4 0.1162 0.1162 1 181.096
## age:location:salinity 5 0.0155 0.0155 1 184.257
## age:location:elevation_sc 6 0.0170 0.0170 1 182.359
## date_planted_grp 7 0.2300 0.2300 1 190.751
## location:co2:salinity:elevation_sc 8 0.5739 0.5739 1 201.997
## location:co2:elevation_sc 9 0.0440 0.0440 1 192.385
## location:co2:salinity 10 0.2111 0.2111 1 184.830
## location:salinity:elevation_sc 11 0.3550 0.3550 1 189.972
## location:salinity 12 0.0194 0.0194 1 192.593
## location:elevation_sc 13 0.1527 0.1527 1 189.981
## age:location:co2 14 0.6496 0.6496 1 196.121
## location:co2 15 0.2865 0.2865 1 196.906
## age:location 16 0.5040 0.5040 1 24.407
## weight_init_sc 0 1.7330 1.7330 1 207.800
## location 0 0.0558 0.0558 1 25.858
## I(elevation_sc^2) 0 6.4346 6.4346 1 198.512
## age:co2:salinity:elevation_sc 0 1.8198 1.8198 1 206.044
## F value Pr(>F)
## origin_lab 0.0386 0.845437
## age:location:co2:elevation_sc 0.2781 0.598571
## age:location:co2:salinity 0.3622 0.548069
## age:location:salinity:elevation_sc 0.4412 0.507405
## age:location:salinity 0.0588 0.808599
## age:location:elevation_sc 0.0651 0.798966
## date_planted_grp 0.8840 0.348286
## location:co2:salinity:elevation_sc 2.2112 0.138573
## location:co2:elevation_sc 0.1666 0.683634
## location:co2:salinity 0.8054 0.370644
## location:salinity:elevation_sc 1.3544 0.245970
## location:salinity 0.0740 0.785875
## location:elevation_sc 0.5838 0.445767
## age:location:co2 2.4948 0.115837
## location:co2 1.0872 0.298363
## age:location 1.9062 0.179896
## weight_init_sc 6.5577 0.011153 *
## location 0.2110 0.649789
## I(elevation_sc^2) 24.3478 1.698e-06 ***
## age:co2:salinity:elevation_sc 6.8859 0.009338 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## sqrt(agb_scam) ~ weight_init_sc + age + location + co2 + salinity +
## elevation_sc + I(elevation_sc^2) + (1 | genotype) + age:co2 +
## age:salinity + age:elevation_sc + co2:salinity + co2:elevation_sc +
## salinity:elevation_sc + age:co2:salinity + age:co2:elevation_sc +
## age:salinity:elevation_sc + co2:salinity:elevation_sc + age:co2:salinity:elevation_sc
# Fit the proposed stepwise model
agb_mod4 <- lmer(sqrt(agb_scam) ~ weight_init + I(elevation^2) + location +
age*co2*salinity*elevation + (1 | genotype), data = bg_nocomp)
# See if there are random effects models that are an improvement of this model
agb_mod5 <- update(agb_mod4, .~.-(1|genotype) + (co2*elevation_sc|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod6 <- update(agb_mod4, .~.-(1|genotype) + (co2*salinity|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod7 <- update(agb_mod4, .~.-(1|genotype) + (salinity*elevation_sc|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod8 <- update(agb_mod4, .~.-(1|genotype) + (co2+elevation_sc|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod9 <- update(agb_mod4, .~.-(1|genotype) + (co2+salinity|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod10 <- update(agb_mod4, .~.-(1|genotype) + (salinity+elevation_sc|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod11 <- update(agb_mod4, .~.-(1|genotype) + (co2|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod12 <- update(agb_mod4, .~.-(1|genotype) + (salinity|genotype))
## boundary (singular) fit: see ?isSingular
agb_mod13 <- update(agb_mod4, .~.-(1|genotype) + (elevation_sc|genotype))
## boundary (singular) fit: see ?isSingular
# So final model is the same as the stepwise proposed model
agb_final <- agb_mod4
# Look at model diagnostics
plot_model(agb_final, type = "diag")
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

##
## [[2]]
## [[2]]$genotype
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[3]]

##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

# Look at random effects
plot_model(agb_final, type = "re")

# Refit model and get Type 3 Anova table
Anova(agb_final, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: sqrt(agb_scam)
## Chisq Df Pr(>Chisq)
## (Intercept) 9.1719 1 0.0024576 **
## weight_init 6.5577 1 0.0104433 *
## I(elevation^2) 24.3478 1 8.042e-07 ***
## location 0.2110 1 0.6459501
## age 2.0667 1 0.1505426
## co2 0.0603 1 0.8059773
## salinity 2.4430 1 0.1180472
## elevation 11.5307 1 0.0006845 ***
## age:co2 5.8183 1 0.0158601 *
## age:salinity 3.7727 1 0.0520955 .
## co2:salinity 3.5302 1 0.0602594 .
## age:elevation 1.8766 1 0.1707246
## co2:elevation 0.1410 1 0.7072758
## salinity:elevation 0.9067 1 0.3409938
## age:co2:salinity 9.8915 1 0.0016604 **
## age:co2:elevation 4.3381 1 0.0372683 *
## age:salinity:elevation 2.2560 1 0.1330943
## co2:salinity:elevation 2.9289 1 0.0870051 .
## age:co2:salinity:elevation 6.8859 1 0.0086878 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Extract predicted means for model
agb_means <- emmeans::emmeans(agb_final, ~co2:salinity:age:elevation, at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")
agb_means_out <- summary(agb_means)
agb_means_out %>%
mutate(salinity = case_when(salinity == "fresh" ~ "Freshwater Site (6ppt)",
T ~ "Brackish Site (8ppt)"),
age = case_when(age == "modern" ~ "Descendant cohort (2000-2020)",
T ~ "Ancestral cohort (1900-1960)")) -> agb_means_out
bg_nocomp_plot <- bg_nocomp %>%
mutate(salinity = case_when(salinity == "fresh" ~ "Freshwater Site (6ppt)",
T ~ "Brackish Site (8ppt)"),
age = case_when(age == "modern" ~ "Descendant cohort (2000-2020)",
T ~ "Ancestral cohort (1900-1960)"))
agb_means_out %>%
ggplot(aes(x = elevation, y = response, color = co2)) +
geom_line(size = 1.5) +
#geom_line(aes(x = elevation, y = lower.CL), linetype = "dashed") +
#geom_line(aes(x = elevation, y = upper.CL), linetype = "dashed") +
facet_grid(age~salinity) +
geom_point(data = bg_nocomp_plot, aes(x = elevation, y = agb_scam), alpha = 0.5) +
theme_bw() +
scale_color_manual(values = c("gray47", "#66a61e")) +
ylab("Aboveground biomass (g)") +
xlab("Elevation (m NAVD88)") +
guides(color=guide_legend(title=expression(CO[2]))) +
theme(legend.position = "top")-> a
# Comparison to a null model that does not include age, location or genotype
agb_null <- lm(sqrt(agb_scam) ~ weight_init_sc +
(co2 + salinity + elevation_sc)^3 + I(elevation_sc^2), data = bg_nocomp)
# Get R2 from model
agb_null_r2 <- r.squaredGLMM(agb_null)[2]
## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.
# Compare to conditional R2 from full model
agb_r2 <- r.squaredGLMM(agb_final)[2]
# Create enough colors for each genotype to have one
n <- length(unique(bg_nocomp$genotype))
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))
# Create caterpillar graph
REsim(agb_final) %>%
mutate(groupID = fct_reorder(groupID, mean)) %>%
ggplot(aes(x = groupID)) +
geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
coord_flip() +
ylab("Random effect") +
xlab("Genotype") +
geom_hline(aes(yintercept = 0), linetype = "dashed") +
theme_bw() + scale_x_discrete(position = "top") +
scale_color_manual(values = col_vector) +
theme(legend.position = "none")-> b
#theme(axis.text.y=element_blank(),
#axis.ticks.y = element_blank()) -> b
cowplot::plot_grid(a,b, rel_widths = c(4,1)) -> agb_plot
Root:shoot ratio
# Calculate root-to-shoot ratio on no comp pots
bg_nocomp %>%
mutate(rs = total_bg / agb_scam) -> bg_nocomp
# Plot and see if there are outliers
bg_nocomp %>%
ggplot(aes(x = rs)) +
geom_histogram() +
geom_rug() +
xlab("root-to-shoot ratio")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite values (stat_bin).

# There are two pots that have really high r:s and one that seems a bit high.
# Take a look at these further.
bg_nocomp %>%
filter(rs > 4.5)
## pot_no beta bg_scam10 bg_sppa10 total_bg site_frame level position
## 1 98 0.9392815 1.460 NA 5.552 fresh_3 1 2
## 2 314 NA 0.084 NA 0.084 fresh_7 4 2
## 3 1936 0.9584345 0.151 NA 0.371 gcrew_6 4 6
## co2 salinity comp elevation genotype location depth_seed cohort
## 1 ambient fresh 0 0.5260000 cm8 corn 0 corn_modern
## 2 ambient fresh 0 0.1980000 cm3 corn 0 corn_modern
## 3 ambient salt 0 0.1853333 ca1 corn 20.75 corn_ancestral
## origin_lab date_cloned_grp date_planted_grp weight_init width_init nodes_init
## 1 nd 1 1 0.667 2.29 2
## 2 utk 1 1 0.724 0.70 1
## 3 utk 1 2 0.390 1.44 1
## agb_scam agb_sppa dens_scam_live dens_scam_dead height_sppa height_scam_tot
## 1 1.186 0 12 0 NA 40.08333
## 2 0.007 0 1 0 NA NA
## 3 0.047 0 1 0 NA 35.50000
## height_scam_grn width_scam_mid width_scam_bas prop_twist soil_to_pot
## 1 38.16667 1.583333 1.775 0.4166667 6.0
## 2 NA NA NA NA 2.5
## 3 34.00000 1.900000 2.300 0.0000000 5.0
## id extinct age elevation_sc weight_init_sc rs
## 1 fresh_3_1_2 0 modern 1.392848 -0.7175561 4.681282
## 2 fresh_7_4_2 0 modern -1.274990 -0.6410989 12.000000
## 3 gcrew_6_4_6 0 ancestral -1.378016 -1.0891112 7.893617
# It seems like pot_no 98 is real -- just a lot of bg compared to ag (on level
# 1)
# The other two seem to be more of experimental artifacts (i.e. initial
# propagule weight is driving total belowground weight at the end), so we drop
# those two
bg_nocomp %>%
filter(rs < 6) -> bg_nocomp_rs
# There also appear to be some that are NAs
bg_nocomp %>%
filter(!complete.cases(rs))
## pot_no beta bg_scam10 bg_sppa10 total_bg site_frame level position co2
## 1 205 NA NA NA NA fresh_5 2 5 elevated
## 2 1700 NA NA NA NA gcrew_4 3 4 elevated
## salinity comp elevation genotype location depth_seed cohort
## 1 fresh 0 0.3973333 ca2 corn 18.75 corn_ancestral
## 2 salt 0 0.2871667 ca4 corn 20.75 corn_ancestral
## origin_lab date_cloned_grp date_planted_grp weight_init width_init nodes_init
## 1 utk 1 1 0.733 1.71 1
## 2 utk 1 2 0.656 1.22 2
## agb_scam agb_sppa dens_scam_live dens_scam_dead height_sppa height_scam_tot
## 1 7.287000 0 32 2 NA 47.51562
## 2 8.005503 0 27 3 NA 48.67308
## height_scam_grn width_scam_mid width_scam_bas prop_twist soil_to_pot
## 1 45.42188 3.096875 3.287500 0.4375000 4.5
## 2 47.13462 3.326923 3.896154 0.3461538 4.0
## id extinct age elevation_sc weight_init_sc rs
## 1 fresh_5_2_5 0 ancestral 0.3463181 -0.6290267 NA
## 2 gcrew_4_3_4 0 ancestral -0.5497391 -0.7323110 NA
# Need to figure this out on data checking!!! I think these were pots with
# messed up bg processing. Move forward with this for now, but I'll go back and
# fix these later
# Change names of salinity factor for presentation
bg_nocomp_rs %>%
mutate(salinity = case_when(salinity == "fresh" ~ "Freshwater Site (6ppt)",
T ~ "Brackish Site (8ppt)")) -> bg_nocomp_rs
# Scale elevation
bg_nocomp_rs$elevation_sc <- scale(bg_nocomp_rs$elevation)[,1]
centering <- attr(scale(bg_nocomp_rs$elevation), "scaled:center")
scaling <- attr(scale(bg_nocomp_rs$elevation), "scaled:scale")
## Fit model and model selection (fixed effects) ####
# Most complex model -- log is the best transformation
rs_mod <- lmer(log(rs) ~ weight_init + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation)^5 + I(elevation^2) +
+ (1|site_frame) + (1|genotype), data = bg_nocomp_rs)
# Do stepwise model selection
# step(rs_mod, keep = c("co2", "salinity", "elevation", "age", "location"))
# Error because model is too complex -- Try with just 4-way interactions
# Most complex model -- log is the best transformation
rs_mod2 <- lmer(log(rs) ~ weight_init + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation)^4 + I(elevation^2) +
+ (1|site_frame) + (1|genotype), data = bg_nocomp_rs)
# Do stepwise model selection
# step(rs_mod2, keep = c("co2", "salinity", "elevation", "age", "location")) # Still mad
Anova(rs_mod2) # Could drop down to 3-way interactions
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: log(rs)
## Chisq Df Pr(>Chisq)
## weight_init 0.4123 1 0.5207998
## date_planted_grp 0.9076 1 0.3407395
## origin_lab 0.2156 1 0.6424063
## age 0.0785 1 0.7792791
## location 0.9893 1 0.3199205
## co2 6.5209 1 0.0106614 *
## salinity 1.5615 1 0.2114511
## elevation 17.5362 1 2.819e-05 ***
## I(elevation^2) 3.5973 1 0.0578750 .
## age:location 0.1646 1 0.6849751
## age:co2 0.0388 1 0.8438226
## age:salinity 0.5912 1 0.4419611
## age:elevation 1.7376 1 0.1874398
## location:co2 0.0051 1 0.9432056
## location:salinity 2.9562 1 0.0855463 .
## location:elevation 5.1197 1 0.0236561 *
## co2:salinity 0.5831 1 0.4451132
## co2:elevation 0.0537 1 0.8166909
## salinity:elevation 13.8903 1 0.0001938 ***
## age:location:co2 0.6379 1 0.4244705
## age:location:salinity 1.6866 1 0.1940492
## age:location:elevation 0.5785 1 0.4468985
## age:co2:salinity 0.0715 1 0.7891903
## age:co2:elevation 5.4149 1 0.0199660 *
## age:salinity:elevation 1.2144 1 0.2704696
## location:co2:salinity 0.0891 1 0.7652656
## location:co2:elevation 0.0846 1 0.7711070
## location:salinity:elevation 2.1855 1 0.1393118
## co2:salinity:elevation 0.1411 1 0.7071706
## age:location:co2:salinity 0.1081 1 0.7423733
## age:location:co2:elevation 0.5091 1 0.4755231
## age:location:salinity:elevation 0.0940 1 0.7591596
## age:co2:salinity:elevation 0.0103 1 0.9190293
## location:co2:salinity:elevation 0.0878 1 0.7669450
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rs_mod3 <- lmer(log(rs) ~ weight_init + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation)^3 + I(elevation^2) +
+ (1|site_frame) + (1|genotype), data = bg_nocomp_rs)
step(rs_mod3, keep = c("co2", "salinity", "elevation", "age", "location"))
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 33 -23.116 112.23
## (1 | site_frame) 1 32 -23.367 110.73 0.5025 1 0.4784
## (1 | genotype) 0 31 -31.591 125.18 16.4477 1 5.001e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## location:co2:elevation 1 0.00470 0.00470 1 182.373 0.0883
## location:co2:salinity 2 0.00699 0.00699 1 176.065 0.1322
## co2:salinity:elevation 3 0.00646 0.00646 1 190.483 0.1227
## age:co2:salinity 4 0.00769 0.00769 1 177.362 0.1470
## origin_lab 5 0.01622 0.01622 1 34.870 0.3115
## weight_init 6 0.02904 0.02904 1 194.177 0.5569
## age:location:co2 7 0.03373 0.03373 1 186.721 0.6473
## location:co2 8 0.00016 0.00016 1 187.020 0.0030
## co2:salinity 9 0.03929 0.03929 1 181.555 0.7608
## age:location:elevation 10 0.04495 0.04495 1 186.711 0.8727
## age:salinity:elevation 11 0.05752 0.05752 1 187.427 1.1129
## date_planted_grp 12 0.05223 0.05223 1 193.891 1.0090
## age:location:salinity 13 0.08469 0.08469 1 191.077 1.6389
## age:location 14 0.00834 0.00834 1 25.358 0.1610
## age:salinity 15 0.03898 0.03898 1 193.725 0.7523
## location:salinity:elevation 16 0.09355 0.09355 1 193.763 1.8097
## location:salinity 17 0.11842 0.11842 1 193.828 2.2880
## I(elevation^2) 18 0.17629 0.17629 1 196.678 3.3896
## location:elevation 0 0.26582 0.26582 1 194.340 5.0620
## salinity:elevation 0 0.71603 0.71603 1 194.686 13.6353
## age:co2:elevation 0 0.28637 0.28637 1 199.725 5.4533
## Pr(>F)
## location:co2:elevation 0.7666359
## location:co2:salinity 0.7165510
## co2:salinity:elevation 0.7265044
## age:co2:salinity 0.7018869
## origin_lab 0.5803156
## weight_init 0.4564266
## age:location:co2 0.4221042
## location:co2 0.9560358
## co2:salinity 0.3842361
## age:location:elevation 0.3514056
## age:salinity:elevation 0.2928045
## date_planted_grp 0.3164039
## age:location:salinity 0.2020318
## age:location 0.6915659
## age:salinity 0.3868097
## location:salinity:elevation 0.1801213
## location:salinity 0.1320070
## I(elevation^2) 0.0671155 .
## location:elevation 0.0255772 *
## salinity:elevation 0.0002882 ***
## age:co2:elevation 0.0205239 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## log(rs) ~ age + location + co2 + salinity + elevation + (1 |
## genotype) + age:co2 + age:elevation + location:elevation +
## co2:elevation + salinity:elevation + age:co2:elevation
# Create the stepwise model
rs_mod4 <- lmer(log(rs) ~ location*elevation + salinity*elevation + age*co2*elevation+
(1|genotype), data = bg_nocomp_rs)
# Diagnostic plots
plot_model(rs_mod4, type = "diag")
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

##
## [[2]]
## [[2]]$genotype
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[3]]

##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

# See if there are random effects models that are an improvement of this model
rs_mod5 <- update(rs_mod4, .~.-(1|genotype) + (co2*elevation|genotype))
## boundary (singular) fit: see ?isSingular
rs_mod6 <- update(rs_mod4, .~.-(1|genotype) + (co2*salinity|genotype))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -8.2e-03
rs_mod7 <- update(rs_mod4, .~.-(1|genotype) + (salinity*elevation|genotype)) # fits
rs_mod8 <- update(rs_mod4, .~.-(1|genotype) + (co2+elevation|genotype))
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -9.4e-05
rs_mod9 <- update(rs_mod4, .~.-(1|genotype) + (co2+salinity|genotype))
rs_mod10 <- update(rs_mod4, .~.-(1|genotype) + (salinity+elevation|genotype)) # fits
rs_mod11 <- update(rs_mod4, .~.-(1|genotype) + (co2|genotype))
## boundary (singular) fit: see ?isSingular
rs_mod12 <- update(rs_mod4, .~.-(1|genotype) + (salinity|genotype)) # fits
rs_mod13 <- update(rs_mod4, .~.-(1|genotype) + (elevation|genotype)) # fits
anova(rs_mod4, rs_mod7, refit = F) # ns
## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age *
## rs_mod4: co2 * elevation + (1 | genotype)
## rs_mod7: log(rs) ~ location + elevation + salinity + age + co2 + (salinity *
## rs_mod7: elevation | genotype) + location:elevation + elevation:salinity +
## rs_mod7: age:co2 + elevation:age + elevation:co2 + elevation:age:co2
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rs_mod4 14 57.793 105.68 -14.897 29.793
## rs_mod7 23 67.283 145.96 -10.642 21.283 8.5103 9 0.4836
anova(rs_mod4, rs_mod10, refit = F) # ns
## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age *
## rs_mod4: co2 * elevation + (1 | genotype)
## rs_mod10: log(rs) ~ location + elevation + salinity + age + co2 + (salinity +
## rs_mod10: elevation | genotype) + location:elevation + elevation:salinity +
## rs_mod10: age:co2 + elevation:age + elevation:co2 + elevation:age:co2
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rs_mod4 14 57.793 105.68 -14.897 29.793
## rs_mod10 19 62.433 127.42 -12.217 24.433 5.36 5 0.3735
anova(rs_mod4, rs_mod12, refit = F) # ns
## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age *
## rs_mod4: co2 * elevation + (1 | genotype)
## rs_mod12: log(rs) ~ location + elevation + salinity + age + co2 + (salinity |
## rs_mod12: genotype) + location:elevation + elevation:salinity + age:co2 +
## rs_mod12: elevation:age + elevation:co2 + elevation:age:co2
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rs_mod4 14 57.793 105.68 -14.897 29.793
## rs_mod12 16 58.109 112.84 -13.055 26.109 3.6842 2 0.1585
anova(rs_mod4, rs_mod13, refit = F) # ns
## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age *
## rs_mod4: co2 * elevation + (1 | genotype)
## rs_mod13: log(rs) ~ location + elevation + salinity + age + co2 + (elevation |
## rs_mod13: genotype) + location:elevation + elevation:salinity + age:co2 +
## rs_mod13: elevation:age + elevation:co2 + elevation:age:co2
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rs_mod4 14 57.793 105.68 -14.897 29.793
## rs_mod13 16 60.732 115.46 -14.366 28.732 1.0617 2 0.5881
# None of the random effects add significantly
# Create caterpillar graph
REsim(rs_mod4) %>%
mutate(groupID = fct_reorder(groupID, mean)) %>%
ggplot(aes(x = groupID)) +
geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
coord_flip() +
ylab("Random effect") +
xlab("Genotype") +
geom_hline(aes(yintercept = 0), linetype = "dashed") +
theme_bw() + scale_x_discrete(position = "top") +
theme(legend.position = "none") +
scale_color_manual(values = col_vector)-> d
# Also fit a null model
rs_mod_null <- lm(log(rs) ~ elevation*salinity + co2*elevation, data = bg_nocomp_rs)
# Compare R2 values
rs_null_r2 <- r.squaredGLMM(rs_mod_null)[2]
rs_r2 <- r.squaredGLMM(rs_mod4)[2]
# Make plots of fixed effects
Anova(rs_mod4, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: log(rs)
## Chisq Df Pr(>Chisq)
## (Intercept) 7.4957 1 0.006185 **
## location 2.1537 1 0.142231
## elevation 56.3481 1 6.071e-14 ***
## salinity 18.6729 1 1.552e-05 ***
## age 0.4679 1 0.493940
## co2 0.4316 1 0.511207
## location:elevation 5.0620 1 0.024456 *
## elevation:salinity 13.6353 1 0.000222 ***
## age:co2 4.9595 1 0.025948 *
## elevation:age 0.3915 1 0.531489
## elevation:co2 2.0213 1 0.155106
## elevation:age:co2 5.4533 1 0.019531 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Need elevation:age:co2, elevation:salinity, and location:elevation
# Extract predicted means for model
rs_means <- emmeans::emmeans(rs_mod4, ~co2:age:elevation,
at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")
rs_means_out <- summary(rs_means)
rs_means_out %>%
mutate(age = case_when(age == "modern" ~ "Descendant cohort (2000-2020)",
T ~ "Ancestral cohort (1900-1960)")) -> rs_means_out
bg_nocomp_rs_plot <- bg_nocomp_plot %>%
mutate(rs = total_bg / agb_scam) %>%
filter(rs < 5)
rs_means_out %>%
ggplot(aes(x = elevation, y = response, color = co2)) +
geom_line(size = 1.5) +
facet_wrap(~age) +
geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = rs), alpha = 0.5) +
theme_bw() +
scale_color_manual(values = c("gray47", "#66a61e")) +
ylab("Root-to-shoot ratio") +
xlab("Elevation (m NAVD88)") +
guides(color=guide_legend(title=expression(CO[2]))) +
theme(legend.position = "top") -> e
rs_means2 <- emmeans(rs_mod4, ~salinity:elevation,
at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")
rs_means_out2 <- summary(rs_means2)
rs_means_out2 %>%
ggplot(aes(x = elevation, y = response, color = salinity)) +
geom_line(size = 1.5) +
geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = rs), alpha = 0.5) +
theme_bw() + scale_color_manual(values = c("orange", "gray47")) +
ylab("Root-to-shoot ratio") +
xlab("Elevation (m NAVD88)") +
guides(color=guide_legend("")) +
theme(legend.position = "top") -> f
rs_means3 <- emmeans(rs_mod4, ~location:elevation,
at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")
rs_means_out3 <- summary(rs_means3)
rs_means_out3 %>%
mutate(location = case_when(location == "corn" ~ "Corn Island",
T ~ "Kirkpatrick Marsh")) -> rs_means_out3
bg_nocomp_rs_plot %>%
mutate(location = case_when(location == "corn" ~ "Corn Island",
T ~ "Kirkpatrick Marsh")) -> bg_nocomp_rs_plot
rs_means_out3 %>%
ggplot(aes(x = elevation, y = response, color = location)) +
geom_line(size = 1.5) +
geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = rs), alpha = 0.5) +
theme_bw() + scale_color_manual(values = c("#1b9e77", "#7570b3")) +
ylab("Root-to-shoot ratio") +
xlab("Elevation (m NAVD88)") +
guides(color=guide_legend("")) +
theme(legend.position = "top") -> g
cowplot::plot_grid(f, g, labels = c("b", "c")) -> fg
cowplot::plot_grid(e, fg, ncol = 1, labels = c("a", "", "")) -> fge
cowplot::plot_grid(d, labels = "d") -> d
fge + d + plot_layout(widths = c(4,1), guides = "collect") -> rs_plot
Root distribution parameter
# Most complex model
beta_mod <- lmer(beta ~ weight_init + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|site_frame) + (1|genotype), data = bg_nocomp_rs)
# Do stepwise model selection
# step(beta_mod) # too complex
VarCorr(beta_mod)
## Groups Name Std.Dev.
## genotype (Intercept) 0.0119866
## site_frame (Intercept) 0.0032349
## Residual 0.0266574
Anova(beta_mod)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: beta
## Chisq Df Pr(>Chisq)
## weight_init 1.8577 1 0.17289
## date_planted_grp 0.2637 1 0.60761
## origin_lab 0.8143 1 0.36686
## age 3.6133 1 0.05732 .
## location 0.3393 1 0.56025
## co2 1.1127 1 0.29150
## salinity 0.5640 1 0.45264
## elevation_sc 204.6761 1 < 2e-16 ***
## I(elevation_sc^2) 2.4456 1 0.11785
## age:location 0.0324 1 0.85722
## age:co2 0.2908 1 0.58973
## age:salinity 1.0714 1 0.30062
## age:elevation_sc 1.5983 1 0.20614
## location:co2 0.0435 1 0.83487
## location:salinity 4.0349 1 0.04457 *
## location:elevation_sc 4.2626 1 0.03896 *
## co2:salinity 0.5169 1 0.47215
## co2:elevation_sc 4.2644 1 0.03892 *
## salinity:elevation_sc 5.9665 1 0.01458 *
## age:location:co2 1.1063 1 0.29288
## age:location:salinity 0.6866 1 0.40731
## age:location:elevation_sc 0.3536 1 0.55207
## age:co2:salinity 0.0112 1 0.91580
## age:co2:elevation_sc 2.8633 1 0.09062 .
## age:salinity:elevation_sc 1.2128 1 0.27079
## location:co2:salinity 0.1042 1 0.74688
## location:co2:elevation_sc 0.4084 1 0.52277
## location:salinity:elevation_sc 1.0313 1 0.30984
## co2:salinity:elevation_sc 2.7875 1 0.09500 .
## age:location:co2:salinity 0.0041 1 0.94911
## age:location:co2:elevation_sc 0.6246 1 0.42933
## age:location:salinity:elevation_sc 1.4932 1 0.22172
## age:co2:salinity:elevation_sc 0.1184 1 0.73074
## location:co2:salinity:elevation_sc 0.1714 1 0.67886
## age:location:co2:salinity:elevation_sc 0.1397 1 0.70859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# site_frame variance is low and no significant 5-way and 4-way interactions
beta_mod2 <- update(beta_mod, .~.-(1|site_frame))
# step(beta_mod2)
beta_mod3 <- lmer(beta ~ weight_init + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^4 + I(elevation_sc^2) +
(1|genotype), data = bg_nocomp_rs)
step(beta_mod3, keep = c("co2", "salinity", "elevation_sc", "age", "location"))
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 37 355.15 -636.30
## (1 | genotype) 0 36 350.49 -628.98 9.3233 1 0.002263 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## age:location:co2:salinity 1 0.0000035 0.0000035 1 169.793
## age:co2:salinity:elevation_sc 2 0.0000897 0.0000897 1 182.482
## age:co2:salinity 3 0.0000147 0.0000147 1 169.454
## location:co2:salinity:elevation_sc 4 0.0001251 0.0001251 1 190.543
## location:co2:salinity 5 0.0000678 0.0000678 1 171.902
## date_planted_grp 6 0.0001786 0.0001786 1 183.350
## origin_lab 7 0.0005246 0.0005246 1 33.553
## age:location:co2:elevation_sc 8 0.0004973 0.0004973 1 182.509
## location:co2:elevation_sc 9 0.0002324 0.0002324 1 182.652
## age:location:co2 10 0.0008218 0.0008218 1 182.753
## location:co2 11 0.0000119 0.0000119 1 183.807
## age:location:salinity:elevation_sc 12 0.0008337 0.0008337 1 183.998
## age:location:elevation_sc 13 0.0001855 0.0001855 1 183.555
## age:location:salinity 14 0.0004579 0.0004579 1 184.733
## age:location 15 0.0000516 0.0000516 1 26.153
## location:salinity:elevation_sc 16 0.0007878 0.0007878 1 188.462
## age:salinity:elevation_sc 17 0.0008435 0.0008435 1 186.440
## age:salinity 18 0.0005978 0.0005978 1 189.759
## weight_init 19 0.0007834 0.0007834 1 201.412
## co2:salinity:elevation_sc 20 0.0013889 0.0013889 1 200.890
## co2:salinity 21 0.0006082 0.0006082 1 186.612
## I(elevation_sc^2) 22 0.0015584 0.0015584 1 194.340
## age:co2:elevation_sc 23 0.0025997 0.0025997 1 197.628
## age:co2 24 0.0002661 0.0002661 1 197.664
## age:elevation_sc 25 0.0013573 0.0013573 1 192.175
## co2:elevation_sc 26 0.0020924 0.0020924 1 199.335
## location:salinity 27 0.0025404 0.0025404 1 196.318
## location:elevation_sc 28 0.0026331 0.0026331 1 195.520
## age 0 0.0026093 0.0026093 1 26.707
## location 0 0.0005226 0.0005226 1 27.041
## co2 0 0.0008293 0.0008293 1 202.669
## salinity:elevation_sc 0 0.0039719 0.0039719 1 198.555
## F value Pr(>F)
## age:location:co2:salinity 0.0049 0.94400
## age:co2:salinity:elevation_sc 0.1263 0.72267
## age:co2:salinity 0.0208 0.88561
## location:co2:salinity:elevation_sc 0.1783 0.67336
## location:co2:salinity 0.0972 0.75558
## date_planted_grp 0.2575 0.61242
## origin_lab 0.7591 0.38981
## age:location:co2:elevation_sc 0.7203 0.39716
## location:co2:elevation_sc 0.3374 0.56204
## age:location:co2 1.1999 0.27478
## location:co2 0.0172 0.89568
## age:location:salinity:elevation_sc 1.2178 0.27123
## age:location:elevation_sc 0.2713 0.60309
## age:location:salinity 0.6732 0.41298
## age:location 0.0761 0.78482
## location:salinity:elevation_sc 1.1625 0.28232
## age:salinity:elevation_sc 1.2465 0.26566
## age:salinity 0.8816 0.34895
## weight_init 1.1559 0.28360
## co2:salinity:elevation_sc 2.0460 0.15416
## co2:salinity 0.8908 0.34648
## I(elevation_sc^2) 2.2878 0.13202
## age:co2:elevation_sc 3.7853 0.05312 .
## age:co2 0.3853 0.53548
## age:elevation_sc 1.9713 0.16192
## co2:elevation_sc 3.0235 0.08361 .
## location:salinity 3.6577 0.05727 .
## location:elevation_sc 3.7608 0.05391 .
## age 3.6662 0.06629 .
## location 0.7343 0.39902
## co2 1.1652 0.28167
## salinity:elevation_sc 5.5807 0.01913 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## beta ~ age + location + co2 + salinity + elevation_sc + (1 |
## genotype) + salinity:elevation_sc
# Create stepwise model
beta_mod4 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 | genotype), data = bg_nocomp_rs)
# Check diagnostics
plot_model(beta_mod4, type = "diag")
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

##
## [[2]]
## [[2]]$genotype
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[3]]

##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

# Some evidence of non-normality, but not horrible
# See if there are significant random intercepts
beta_mod5 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + salinity*elevation_sc| genotype), data = bg_nocomp_rs) # fits
beta_mod6 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + salinity+elevation_sc| genotype), data = bg_nocomp_rs) # fits
beta_mod7 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + salinity+co2| genotype), data = bg_nocomp_rs)
## boundary (singular) fit: see ?isSingular
beta_mod8 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + elevation_sc+co2| genotype), data = bg_nocomp_rs)
## boundary (singular) fit: see ?isSingular
beta_mod9 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + co2| genotype), data = bg_nocomp_rs)
## boundary (singular) fit: see ?isSingular
beta_mod10 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + salinity| genotype), data = bg_nocomp_rs) # fits
beta_mod11 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + elevation_sc| genotype), data = bg_nocomp_rs) # fits
anova(beta_mod5, beta_mod4) # *** best model
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp_rs
## Models:
## beta_mod4: beta ~ age + location + co2 + salinity * elevation_sc + (1 |
## beta_mod4: genotype)
## beta_mod5: beta ~ age + location + co2 + salinity * elevation_sc + (1 +
## beta_mod5: salinity * elevation_sc | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## beta_mod4 9 -940.68 -910.02 479.34 -958.68
## beta_mod5 18 -958.06 -896.73 497.03 -994.06 35.371 9 5.124e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(beta_mod6, beta_mod5) # ***
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp_rs
## Models:
## beta_mod6: beta ~ age + location + co2 + salinity * elevation_sc + (1 +
## beta_mod6: salinity + elevation_sc | genotype)
## beta_mod5: beta ~ age + location + co2 + salinity * elevation_sc + (1 +
## beta_mod5: salinity * elevation_sc | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## beta_mod6 14 -955.93 -908.23 491.97 -983.93
## beta_mod5 18 -958.06 -896.73 497.03 -994.06 10.122 4 0.03843 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(beta_mod10, beta_mod4) # *
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp_rs
## Models:
## beta_mod4: beta ~ age + location + co2 + salinity * elevation_sc + (1 |
## beta_mod4: genotype)
## beta_mod10: beta ~ age + location + co2 + salinity * elevation_sc + (1 +
## beta_mod10: salinity | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## beta_mod4 9 -940.68 -910.02 479.34 -958.68
## beta_mod10 11 -943.92 -906.44 482.96 -965.92 7.2339 2 0.02686 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(beta_mod11, beta_mod4) # **
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp_rs
## Models:
## beta_mod4: beta ~ age + location + co2 + salinity * elevation_sc + (1 |
## beta_mod4: genotype)
## beta_mod11: beta ~ age + location + co2 + salinity * elevation_sc + (1 +
## beta_mod11: elevation_sc | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## beta_mod4 9 -940.68 -910.02 479.34 -958.68
## beta_mod11 11 -946.63 -909.15 484.31 -968.63 9.9407 2 0.006941 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Best model has salinity*elevation random slope
beta_finalmod <-lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + salinity*elevation_sc| genotype), data = bg_nocomp_rs)
# Look at random effect
out <- plot_model(beta_finalmod, type = "pred", pred.type = "re", terms = c("elevation_sc[all]", "genotype","salinity")) +
scale_color_manual(values = rainbow(31)) +
ggtitle("") +
theme_bw() +
theme(legend.position = "none") +
ylab("Root distribution parameter") +
xlab("Elevation (m NAVD88)")
## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Set1 is 9
## Returning the palette you asked for with that many colors
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
# Extract values from plot_model so we can customize
tibble(elevation_sc = out$data$x,
root_param = out$data$predicted,
genotype = out$data$group,
salinity = out$data$facet) %>%
mutate(elevation = elevation_sc*scaling + centering) -> beta_GxE_data
beta_GxE_data %>%
filter(elevation %in% min(elevation)) -> mins
beta_GxE_data %>%
filter(elevation %in% max(elevation)) -> maxs
# Make random effects plot
beta_GxE_data %>%
ggplot(aes(x = elevation, y = root_param, color = genotype)) +
geom_line(size = 1, alpha = 1) +
facet_wrap(~salinity) +
scale_color_manual(values = col_vector) +
theme_bw() +
theme(legend.position = "none") +
ylab("Root distribution parameter") +
xlab("Elevation (m NAVD88)") +
geom_point(data = mins, aes(x = elevation, y = root_param, color = genotype), size = 2.7) +
geom_point(data = maxs, aes(x = elevation, y = root_param, color = genotype), size = 2.7) -> h
# Make fixed effect plot
beta_means <- emmeans(beta_finalmod, ~salinity:elevation_sc,
at= list(elevation_sc = seq(-1.72, 1.54, 0.1)), type = "response")
beta_means_out <- summary(beta_means)
beta_means_out %>%
mutate(elevation = elevation_sc*scaling + centering) %>%
ggplot(aes(x = elevation, y = emmean, color = salinity)) +
geom_line(size = 1.5) +
geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = beta), alpha = 0.5) +
theme_bw() + scale_color_manual(values = c("orange", "gray47")) +
ylab("Root distribution parameter") +
xlab("Elevation (m NAVD88)") +
guides(color=guide_legend("")) +
theme(legend.position = "top") +
annotate("text", x = c(0.22, 0.22), y = c(0.95, 0.77), label = c("deeper rooting", "shallower rooting"),color = "gray47", fontface = 3) -> i
beta_finalmod <-lmer(beta ~ age + location + co2 + salinity*elevation_sc +
(1 + salinity*elevation_sc| genotype), data = bg_nocomp_rs_plot)
my_labeller <- as_labeller(function(x){
return(paste0(x))
})
# Make fixed effect plot for age overall
emmip(beta_finalmod, location ~ age, CI = T, weights = "proportional") +
geom_jitter(data = bg_nocomp_rs_plot, aes(x = age, y = beta, color = location),
height = 0, width = 0.1, alpha = 0.3, shape = 1) +
ylab("Root distribution parameter") +
scale_x_discrete("Age cohort",
labels = c("Descendant cohort (2000-2020)" = "Descendant (2000-2020)",
"Ancestral cohort (1900-1960)" = "Ancestral (1900-1960)")) +
scale_color_manual(values = c("#1b9e77", "#7570b3")) + theme_bw() +
theme(legend.position = "top") + guides(color=guide_legend("")) +
scale_size_manual(values=c(2,2)) -> j
(i + j) / h + plot_annotation(tag_levels = 'a') -> beta_plot
# Compare to a null model
beta_null_mod <-lm(beta ~ co2 + salinity*elevation_sc, data = bg_nocomp_rs_plot)
beta_null_r2 <- r.squaredGLMM(beta_null_mod)[2]
beta_r2 <- r.squaredGLMM(beta_finalmod)[2]
Average stem height
height_mod <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|site_frame) + (1|genotype), data = bg_nocomp)
plot(height_mod)# Major outlier

bg_nocomp %>% filter(pot_no == 1830)
## pot_no beta bg_scam10 bg_sppa10 total_bg site_frame level position
## 1 1830 0.93414 0.637 NA 1.239 gcrew_5 1 2
## co2 salinity comp elevation genotype location depth_seed cohort
## 1 elevated salt 0 0.5263333 ca7 corn 20.75 corn_ancestral
## origin_lab date_cloned_grp date_planted_grp weight_init width_init nodes_init
## 1 utk 1 2 0.948 1.69 1
## agb_scam agb_sppa dens_scam_live dens_scam_dead height_sppa height_scam_tot
## 1 0.721 0 7 1 NA 11.5
## height_scam_grn width_scam_mid width_scam_bas prop_twist soil_to_pot
## 1 11.5 1.7 1.7 0 6
## id extinct age elevation_sc weight_init_sc rs
## 1 gcrew_5_1_2 0 ancestral 1.395559 -0.3406356 1.718447
# This is the outlier and note on data sheet says that this pot had significant
# herbivory which would affect the average height
bg_nocomp_height <- bg_nocomp %>%
filter(pot_no != 1830)
height_mod2 <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|site_frame) + (1|genotype), data = bg_nocomp_height)
# Stepwise selection
#step(height_mod2) # Too complex
VarCorr(height_mod2)
## Groups Name Std.Dev.
## genotype (Intercept) 3.57274
## site_frame (Intercept) 0.34105
## Residual 4.62926
Anova(height_mod2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: height_scam_tot
## Chisq Df Pr(>Chisq)
## weight_init_sc 2.5785 1 0.10832
## date_planted_grp 2.5936 1 0.10730
## origin_lab 0.1670 1 0.68278
## age 0.1072 1 0.74338
## location 0.0539 1 0.81645
## co2 0.0164 1 0.89818
## salinity 0.1062 1 0.74455
## elevation_sc 147.9034 1 < 2e-16 ***
## I(elevation_sc^2) 0.2943 1 0.58747
## age:location 0.5482 1 0.45904
## age:co2 0.5388 1 0.46291
## age:salinity 0.5105 1 0.47494
## age:elevation_sc 0.9837 1 0.32128
## location:co2 0.7792 1 0.37739
## location:salinity 3.8988 1 0.04832 *
## location:elevation_sc 2.4535 1 0.11726
## co2:salinity 0.3058 1 0.58028
## co2:elevation_sc 0.8110 1 0.36781
## salinity:elevation_sc 0.5857 1 0.44408
## age:location:co2 0.8210 1 0.36488
## age:location:salinity 0.1151 1 0.73437
## age:location:elevation_sc 1.4778 1 0.22412
## age:co2:salinity 0.1942 1 0.65940
## age:co2:elevation_sc 1.0207 1 0.31235
## age:salinity:elevation_sc 1.1686 1 0.27968
## location:co2:salinity 0.3695 1 0.54330
## location:co2:elevation_sc 0.2627 1 0.60830
## location:salinity:elevation_sc 0.1723 1 0.67808
## co2:salinity:elevation_sc 0.7665 1 0.38131
## age:location:co2:salinity 1.3312 1 0.24859
## age:location:co2:elevation_sc 1.7941 1 0.18043
## age:location:salinity:elevation_sc 0.4387 1 0.50776
## age:co2:salinity:elevation_sc 1.1499 1 0.28358
## location:co2:salinity:elevation_sc 0.3212 1 0.57086
## age:location:co2:salinity:elevation_sc 0.5296 1 0.46677
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Drop 5-way interaction
height_mod3 <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^4 + I(elevation_sc^2) +
(1|site_frame) + (1|genotype), data = bg_nocomp_height)
#step(height_mod3)
VarCorr(height_mod3)
## Groups Name Std.Dev.
## genotype (Intercept) 3.60979
## site_frame (Intercept) 0.35074
## Residual 4.61663
height_mod4 <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^4 +
+I(elevation_sc^2) + (1|genotype), data = bg_nocomp_height)
step(height_mod4, keep = c("co2", "salinity", "age", "location", "elevation_sc"))
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 37 -638.59 1351.2
## (1 | genotype) 0 36 -662.91 1397.8 48.642 1 3.072e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## origin_lab 1 3.4 3.4 1 42.738
## I(elevation_sc^2) 2 5.6 5.6 1 177.115
## location:co2:salinity:elevation_sc 3 6.9 6.9 1 189.353
## age:location:salinity:elevation_sc 4 7.2 7.2 1 176.515
## location:salinity:elevation_sc 5 3.4 3.4 1 177.694
## age:co2:salinity:elevation_sc 6 24.0 24.0 1 186.491
## co2:salinity:elevation_sc 7 14.2 14.2 1 186.332
## age:salinity:elevation_sc 8 22.0 22.0 1 179.905
## salinity:elevation_sc 9 11.2 11.2 1 181.159
## age:location:co2:salinity 10 25.6 25.6 1 180.483
## age:location:salinity 11 2.2 2.2 1 183.906
## age:co2:salinity 12 6.1 6.1 1 182.057
## location:co2:salinity 13 9.9 9.9 1 183.148
## co2:salinity 14 5.8 5.8 1 183.800
## age:salinity 15 12.1 12.1 1 188.737
## age:location:co2:elevation_sc 16 26.5 26.5 1 190.536
## location:co2:elevation_sc 17 4.1 4.1 1 192.002
## age:location:co2 18 17.1 17.1 1 191.656
## location:co2 19 19.5 19.5 1 193.314
## age:co2:elevation_sc 20 25.1 25.1 1 195.232
## age:co2 21 8.4 8.4 1 196.154
## co2:elevation_sc 22 17.6 17.6 1 197.555
## age:location:elevation_sc 23 38.3 38.3 1 194.616
## age:location 24 10.7 10.7 1 28.090
## age:elevation_sc 25 19.2 19.2 1 195.608
## location:elevation_sc 26 31.7 31.7 1 196.096
## weight_init_sc 27 40.1 40.1 1 206.394
## location:salinity 28 70.2 70.2 1 200.231
## date_planted_grp 0 145.9 145.9 1 202.310
## age 0 1.6 1.6 1 29.057
## location 0 2.1 2.1 1 29.310
## co2 0 2.1 2.1 1 203.494
## salinity 0 18.6 18.6 1 202.038
## elevation_sc 0 3429.1 3429.1 1 198.957
## F value Pr(>F)
## origin_lab 0.1588 0.692275
## I(elevation_sc^2) 0.2616 0.609671
## location:co2:salinity:elevation_sc 0.3230 0.570484
## age:location:salinity:elevation_sc 0.3405 0.560301
## location:salinity:elevation_sc 0.1619 0.687876
## age:co2:salinity:elevation_sc 1.1398 0.287079
## co2:salinity:elevation_sc 0.6754 0.412225
## age:salinity:elevation_sc 1.0475 0.307451
## salinity:elevation_sc 0.5337 0.466014
## age:location:co2:salinity 1.2224 0.270353
## age:location:salinity 0.1063 0.744802
## age:co2:salinity 0.2907 0.590431
## location:co2:salinity 0.4762 0.491030
## co2:salinity 0.2781 0.598563
## age:salinity 0.5844 0.445557
## age:location:co2:elevation_sc 1.2846 0.258471
## location:co2:elevation_sc 0.1965 0.658048
## age:location:co2 0.8318 0.362889
## location:co2 0.9496 0.331046
## age:co2:elevation_sc 1.2178 0.271159
## age:co2 0.4051 0.525202
## co2:elevation_sc 0.8533 0.356736
## age:location:elevation_sc 1.8618 0.173992
## age:location 0.5202 0.476699
## age:elevation_sc 0.9287 0.336387
## location:elevation_sc 1.5365 0.216616
## weight_init_sc 1.9409 0.165069
## location:salinity 3.3881 0.067148 .
## date_planted_grp 6.9633 0.008967 **
## age 0.0758 0.784991
## location 0.0986 0.755695
## co2 0.1012 0.750707
## salinity 0.8867 0.347500
## elevation_sc 163.6906 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## height_scam_tot ~ date_planted_grp + age + location + co2 + salinity +
## elevation_sc + (1 | genotype)
# Fit proposed model from step procedure
height_mod5 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
(1 | genotype), data = bg_nocomp_height)
# Try fitting additional random slopes
height_mod6 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
(1+salinity+co2|genotype), data = bg_nocomp_height)
## boundary (singular) fit: see ?isSingular
height_mod7 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
(1+salinity+elevation_sc|genotype), data = bg_nocomp_height)
## boundary (singular) fit: see ?isSingular
height_mod8 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
(1+elevation_sc+co2|genotype), data = bg_nocomp_height) # fit
height_mod9 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
(1+co2|genotype), data = bg_nocomp_height) # fit
height_mod10 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
(1+elevation_sc|genotype), data = bg_nocomp_height) # fit
height_mod11 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
(1+salinity|genotype), data = bg_nocomp_height)
## boundary (singular) fit: see ?isSingular
anova(height_mod5, height_mod8)
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp_height
## Models:
## height_mod5: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity +
## height_mod5: elevation_sc + (1 | genotype)
## height_mod8: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity +
## height_mod8: elevation_sc + (1 + elevation_sc + co2 | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## height_mod5 9 1402.6 1433.4 -692.29 1384.6
## height_mod8 14 1406.2 1454.2 -689.12 1378.2 6.3423 5 0.2743
anova(height_mod5, height_mod9) ## This is the best, but not significant
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp_height
## Models:
## height_mod5: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity +
## height_mod5: elevation_sc + (1 | genotype)
## height_mod9: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity +
## height_mod9: elevation_sc + (1 + co2 | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## height_mod5 9 1402.6 1433.4 -692.29 1384.6
## height_mod9 11 1402.0 1439.8 -690.02 1380.0 4.5389 2 0.1034
anova(height_mod5, height_mod10)
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp_height
## Models:
## height_mod5: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity +
## height_mod5: elevation_sc + (1 | genotype)
## height_mod10: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity +
## height_mod10: elevation_sc + (1 + elevation_sc | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## height_mod5 9 1402.6 1433.4 -692.29 1384.6
## height_mod10 11 1405.7 1443.5 -691.87 1383.7 0.836 2 0.6584
# Create plot of fixed effects
Anova(height_mod5) # Just elevation really
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: height_scam_tot
## Chisq Df Pr(>Chisq)
## date_planted_grp 6.9633 1 0.00832 **
## age 0.0758 1 0.78304
## location 0.0986 1 0.75347
## co2 0.1012 1 0.75038
## salinity 0.8867 1 0.34638
## elevation_sc 163.6906 1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Refit so it is not scaled for easy plotting
height_modfinal <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation +
(1|genotype), data = bg_nocomp_height)
height_means <- emmeans(height_modfinal, ~elevation, at = list(elevation = seq(min(bg_nocomp_height$elevation), max(bg_nocomp_height$elevation), 0.01)), weights = "proportional")
summary(height_means) %>%
ggplot(aes(x = elevation, y = emmean)) +
geom_line(size = 2) +
geom_line(aes(x = elevation, y = lower.CL), linetype = "dashed")+
geom_line(aes(x = elevation, y = upper.CL), linetype = "dashed")+
geom_point(data = bg_nocomp_height, aes(x = elevation, y = height_scam_tot), alpha = 0.2) +
ylab("Mean stem height (cm)") +
xlab("Elevation (m NAVD88)") -> k
# Random effects
REsim(height_modfinal) %>%
mutate(groupID = fct_reorder(groupID, mean)) %>%
ggplot(aes(x = groupID)) +
geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
coord_flip() +
ylab("Random effect") +
xlab("Genotype") +
geom_hline(aes(yintercept = 0), linetype = "dashed") +
theme_bw() + scale_x_discrete(position = "top") +
theme(legend.position = "none") +
scale_color_manual(values = col_vector) -> l
k + l + plot_layout(widths = c(4,1)) +plot_annotation(tag_levels = 'a') -> height_plot
# Fit null model (no age, location or genotype) and compare to height model
height_mod_null <- lm(height_scam_tot ~ date_planted_grp + co2 + salinity + elevation, data = bg_nocomp_height)
height_null_r2 <- r.squaredGLMM(height_mod_null)[2]
height_r2 <- r.squaredGLMM(height_modfinal)[2]
Average stem width
width_mod <- lmer(width_scam_mid ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|site_frame) + (1|genotype), data = bg_nocomp)
plot(width_mod)

#step(width_mod, keep = c("co2", "salinity", "elevation_sc", "age", "location"))
summary(width_mod) # site_frame explains nothing
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_scam_mid ~ weight_init_sc + date_planted_grp + origin_lab +
## (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
## (1 | site_frame) + (1 | genotype)
## Data: bg_nocomp
##
## REML criterion at convergence: 362.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6957 -0.5571 -0.0710 0.5351 2.5515
##
## Random effects:
## Groups Name Variance Std.Dev.
## genotype (Intercept) 0.0909431 0.30157
## site_frame (Intercept) 0.0009338 0.03056
## Residual 0.1914315 0.43753
## Number of obs: 229, groups: genotype, 31; site_frame, 14
##
## Fixed effects:
## Estimate
## (Intercept) 3.552199
## weight_init_sc 0.067583
## date_planted_grp2 -0.405315
## origin_labutk 0.160306
## agemodern -0.236577
## locationkirkpatrick -0.433249
## co2elevated -0.134332
## salinitysalt -0.151769
## elevation_sc -0.469996
## I(elevation_sc^2) 0.006502
## agemodern:locationkirkpatrick 0.449555
## agemodern:co2elevated 0.310725
## agemodern:salinitysalt 0.338239
## agemodern:elevation_sc -0.072590
## locationkirkpatrick:co2elevated 0.206208
## locationkirkpatrick:salinitysalt 0.476190
## locationkirkpatrick:elevation_sc 0.234863
## co2elevated:salinitysalt -0.042617
## co2elevated:elevation_sc 0.045492
## salinitysalt:elevation_sc 0.122281
## agemodern:locationkirkpatrick:co2elevated -0.122541
## agemodern:locationkirkpatrick:salinitysalt -0.369282
## agemodern:locationkirkpatrick:elevation_sc 0.163139
## agemodern:co2elevated:salinitysalt -0.345161
## agemodern:co2elevated:elevation_sc -0.200902
## agemodern:salinitysalt:elevation_sc 0.022314
## locationkirkpatrick:co2elevated:salinitysalt -0.067463
## locationkirkpatrick:co2elevated:elevation_sc -0.248724
## locationkirkpatrick:salinitysalt:elevation_sc -0.189380
## co2elevated:salinitysalt:elevation_sc -0.283756
## agemodern:locationkirkpatrick:co2elevated:salinitysalt 0.134748
## agemodern:locationkirkpatrick:co2elevated:elevation_sc -0.091339
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc -0.023807
## agemodern:co2elevated:salinitysalt:elevation_sc 0.445879
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 0.616211
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc -0.262417
## Std. Error
## (Intercept) 0.180783
## weight_init_sc 0.034249
## date_planted_grp2 0.199876
## origin_labutk 0.134072
## agemodern 0.194852
## locationkirkpatrick 0.234510
## co2elevated 0.170418
## salinitysalt 0.244546
## elevation_sc 0.107438
## I(elevation_sc^2) 0.039324
## agemodern:locationkirkpatrick 0.316275
## agemodern:co2elevated 0.223402
## agemodern:salinitysalt 0.194370
## agemodern:elevation_sc 0.146213
## locationkirkpatrick:co2elevated 0.261749
## locationkirkpatrick:salinitysalt 0.255509
## locationkirkpatrick:elevation_sc 0.182223
## co2elevated:salinitysalt 0.225143
## co2elevated:elevation_sc 0.176680
## salinitysalt:elevation_sc 0.141817
## agemodern:locationkirkpatrick:co2elevated 0.350795
## agemodern:locationkirkpatrick:salinitysalt 0.336443
## agemodern:locationkirkpatrick:elevation_sc 0.251203
## agemodern:co2elevated:salinitysalt 0.304054
## agemodern:co2elevated:elevation_sc 0.241305
## agemodern:salinitysalt:elevation_sc 0.197682
## locationkirkpatrick:co2elevated:salinitysalt 0.379913
## locationkirkpatrick:co2elevated:elevation_sc 0.285121
## locationkirkpatrick:salinitysalt:elevation_sc 0.266749
## co2elevated:salinitysalt:elevation_sc 0.235755
## agemodern:locationkirkpatrick:co2elevated:salinitysalt 0.514829
## agemodern:locationkirkpatrick:co2elevated:elevation_sc 0.408144
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc 0.356388
## agemodern:co2elevated:salinitysalt:elevation_sc 0.324635
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 0.414821
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 0.591532
## df
## (Intercept) 50.726015
## weight_init_sc 157.953703
## date_planted_grp2 179.784008
## origin_labutk 39.580508
## agemodern 54.052246
## locationkirkpatrick 60.502178
## co2elevated 90.361360
## salinitysalt 159.615974
## elevation_sc 168.739610
## I(elevation_sc^2) 171.879127
## agemodern:locationkirkpatrick 56.220115
## agemodern:co2elevated 174.755364
## agemodern:salinitysalt 167.721595
## agemodern:elevation_sc 167.284613
## locationkirkpatrick:co2elevated 170.293590
## locationkirkpatrick:salinitysalt 178.521110
## locationkirkpatrick:elevation_sc 174.608314
## co2elevated:salinitysalt 76.975497
## co2elevated:elevation_sc 168.840090
## salinitysalt:elevation_sc 167.422935
## agemodern:locationkirkpatrick:co2elevated 169.374647
## agemodern:locationkirkpatrick:salinitysalt 176.998758
## agemodern:locationkirkpatrick:elevation_sc 179.011320
## agemodern:co2elevated:salinitysalt 164.962255
## agemodern:co2elevated:elevation_sc 167.803333
## agemodern:salinitysalt:elevation_sc 164.677152
## locationkirkpatrick:co2elevated:salinitysalt 172.679836
## locationkirkpatrick:co2elevated:elevation_sc 177.411207
## locationkirkpatrick:salinitysalt:elevation_sc 177.570616
## co2elevated:salinitysalt:elevation_sc 168.122265
## agemodern:locationkirkpatrick:co2elevated:salinitysalt 169.774274
## agemodern:locationkirkpatrick:co2elevated:elevation_sc 181.543356
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc 179.577296
## agemodern:co2elevated:salinitysalt:elevation_sc 167.919283
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 182.730762
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 184.191479
## t value
## (Intercept) 19.649
## weight_init_sc 1.973
## date_planted_grp2 -2.028
## origin_labutk 1.196
## agemodern -1.214
## locationkirkpatrick -1.847
## co2elevated -0.788
## salinitysalt -0.621
## elevation_sc -4.375
## I(elevation_sc^2) 0.165
## agemodern:locationkirkpatrick 1.421
## agemodern:co2elevated 1.391
## agemodern:salinitysalt 1.740
## agemodern:elevation_sc -0.496
## locationkirkpatrick:co2elevated 0.788
## locationkirkpatrick:salinitysalt 1.864
## locationkirkpatrick:elevation_sc 1.289
## co2elevated:salinitysalt -0.189
## co2elevated:elevation_sc 0.257
## salinitysalt:elevation_sc 0.862
## agemodern:locationkirkpatrick:co2elevated -0.349
## agemodern:locationkirkpatrick:salinitysalt -1.098
## agemodern:locationkirkpatrick:elevation_sc 0.649
## agemodern:co2elevated:salinitysalt -1.135
## agemodern:co2elevated:elevation_sc -0.833
## agemodern:salinitysalt:elevation_sc 0.113
## locationkirkpatrick:co2elevated:salinitysalt -0.178
## locationkirkpatrick:co2elevated:elevation_sc -0.872
## locationkirkpatrick:salinitysalt:elevation_sc -0.710
## co2elevated:salinitysalt:elevation_sc -1.204
## agemodern:locationkirkpatrick:co2elevated:salinitysalt 0.262
## agemodern:locationkirkpatrick:co2elevated:elevation_sc -0.224
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc -0.067
## agemodern:co2elevated:salinitysalt:elevation_sc 1.373
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 1.485
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc -0.444
## Pr(>|t|)
## (Intercept) < 2e-16
## weight_init_sc 0.0502
## date_planted_grp2 0.0441
## origin_labutk 0.2389
## agemodern 0.2300
## locationkirkpatrick 0.0696
## co2elevated 0.4326
## salinitysalt 0.5357
## elevation_sc 2.12e-05
## I(elevation_sc^2) 0.8689
## agemodern:locationkirkpatrick 0.1607
## agemodern:co2elevated 0.1660
## agemodern:salinitysalt 0.0837
## agemodern:elevation_sc 0.6202
## locationkirkpatrick:co2elevated 0.4319
## locationkirkpatrick:salinitysalt 0.0640
## locationkirkpatrick:elevation_sc 0.1991
## co2elevated:salinitysalt 0.8504
## co2elevated:elevation_sc 0.7971
## salinitysalt:elevation_sc 0.3898
## agemodern:locationkirkpatrick:co2elevated 0.7273
## agemodern:locationkirkpatrick:salinitysalt 0.2739
## agemodern:locationkirkpatrick:elevation_sc 0.5169
## agemodern:co2elevated:salinitysalt 0.2579
## agemodern:co2elevated:elevation_sc 0.4063
## agemodern:salinitysalt:elevation_sc 0.9103
## locationkirkpatrick:co2elevated:salinitysalt 0.8593
## locationkirkpatrick:co2elevated:elevation_sc 0.3842
## locationkirkpatrick:salinitysalt:elevation_sc 0.4787
## co2elevated:salinitysalt:elevation_sc 0.2304
## agemodern:locationkirkpatrick:co2elevated:salinitysalt 0.7938
## agemodern:locationkirkpatrick:co2elevated:elevation_sc 0.8232
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc 0.9468
## agemodern:co2elevated:salinitysalt:elevation_sc 0.1714
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 0.1391
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 0.6578
##
## (Intercept) ***
## weight_init_sc .
## date_planted_grp2 *
## origin_labutk
## agemodern
## locationkirkpatrick .
## co2elevated
## salinitysalt
## elevation_sc ***
## I(elevation_sc^2)
## agemodern:locationkirkpatrick
## agemodern:co2elevated
## agemodern:salinitysalt .
## agemodern:elevation_sc
## locationkirkpatrick:co2elevated
## locationkirkpatrick:salinitysalt .
## locationkirkpatrick:elevation_sc
## co2elevated:salinitysalt
## co2elevated:elevation_sc
## salinitysalt:elevation_sc
## agemodern:locationkirkpatrick:co2elevated
## agemodern:locationkirkpatrick:salinitysalt
## agemodern:locationkirkpatrick:elevation_sc
## agemodern:co2elevated:salinitysalt
## agemodern:co2elevated:elevation_sc
## agemodern:salinitysalt:elevation_sc
## locationkirkpatrick:co2elevated:salinitysalt
## locationkirkpatrick:co2elevated:elevation_sc
## locationkirkpatrick:salinitysalt:elevation_sc
## co2elevated:salinitysalt:elevation_sc
## agemodern:locationkirkpatrick:co2elevated:salinitysalt
## agemodern:locationkirkpatrick:co2elevated:elevation_sc
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc
## agemodern:co2elevated:salinitysalt:elevation_sc
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 36 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
width_mod2 <- lmer(width_scam_mid ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|genotype), data = bg_nocomp)
step(width_mod2, keep = c("co2", "salinity", "elevation_sc", "age", "location"))
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 38 -181.14 438.28
## (1 | genotype) 0 37 -198.30 470.60 34.323 1 4.668e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## I(elevation_sc^2) 1 0.00558 0.00558 1 176.642
## age:location:co2:salinity:elevation_sc 2 0.03915 0.03915 1 190.963
## age:location:co2:salinity 3 0.00722 0.00722 1 175.156
## age:location:salinity:elevation_sc 4 0.04394 0.04394 1 178.924
## age:location:co2:elevation_sc 5 0.12646 0.12646 1 184.086
## age:location:co2 6 0.01269 0.01269 1 182.208
## age:location:elevation_sc 7 0.03928 0.03928 1 180.835
## age:location:salinity 8 0.28527 0.28527 1 185.334
## age:location 9 0.21127 0.21127 1 26.253
## origin_lab 10 0.35404 0.35404 1 43.470
## age:co2:salinity:elevation_sc 11 0.38770 0.38770 1 191.499
## age:co2:elevation_sc 12 0.02575 0.02575 1 187.509
## age:salinity:elevation_sc 13 0.22189 0.22189 1 183.576
## age:elevation_sc 14 0.02047 0.02047 1 185.188
## age:co2:salinity 15 0.25857 0.25857 1 183.807
## age:co2 16 0.12120 0.12120 1 191.546
## age:salinity 17 0.16924 0.16924 1 189.549
## location:co2:salinity:elevation_sc 18 0.39122 0.39122 1 205.888
## location:salinity:elevation_sc 19 0.00054 0.00054 1 190.663
## location:co2:elevation_sc 20 0.00836 0.00836 1 193.646
## location:co2:salinity 21 0.01331 0.01331 1 190.110
## co2:salinity:elevation_sc 22 0.08143 0.08143 1 201.648
## location:co2 23 0.14072 0.14072 1 196.758
## co2:elevation_sc 24 0.40056 0.40056 1 197.858
## weight_init_sc 25 0.61674 0.61674 1 204.488
## salinity:elevation_sc 26 0.60474 0.60474 1 197.007
## co2:salinity 27 0.66749 0.66749 1 195.130
## date_planted_grp 0 1.61835 1.61835 1 201.850
## age 0 0.06275 0.06275 1 27.501
## co2 0 0.00276 0.00276 1 203.339
## location:salinity 0 0.81352 0.81352 1 199.950
## location:elevation_sc 0 1.48509 1.48509 1 197.706
## F value Pr(>F)
## I(elevation_sc^2) 0.0290 0.864957
## age:location:co2:salinity:elevation_sc 0.2048 0.651390
## age:location:co2:salinity 0.0379 0.845873
## age:location:salinity:elevation_sc 0.2317 0.630840
## age:location:co2:elevation_sc 0.6701 0.414084
## age:location:co2 0.0675 0.795378
## age:location:elevation_sc 0.2099 0.647398
## age:location:salinity 1.5321 0.217359
## age:location 1.1276 0.297966
## origin_lab 1.8885 0.176419
## age:co2:salinity:elevation_sc 2.0565 0.153194
## age:co2:elevation_sc 0.1361 0.712575
## age:salinity:elevation_sc 1.1785 0.279079
## age:elevation_sc 0.1086 0.742070
## age:co2:salinity 1.3792 0.241763
## age:co2 0.6466 0.422312
## age:salinity 0.9042 0.342882
## location:co2:salinity:elevation_sc 2.0909 0.149695
## location:salinity:elevation_sc 0.0029 0.957164
## location:co2:elevation_sc 0.0447 0.832776
## location:co2:salinity 0.0716 0.789320
## co2:salinity:elevation_sc 0.4402 0.507804
## location:co2 0.7630 0.383471
## co2:elevation_sc 2.1761 0.141757
## weight_init_sc 3.3336 0.069335 .
## salinity:elevation_sc 3.2392 0.073427 .
## co2:salinity 3.5374 0.061490 .
## date_planted_grp 8.4408 0.004078 **
## age 0.3273 0.571898
## co2 0.0144 0.904573
## location:salinity 4.2430 0.040707 *
## location:elevation_sc 7.7457 0.005906 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## width_scam_mid ~ date_planted_grp + age + location + co2 + salinity +
## elevation_sc + (1 | genotype) + location:salinity + location:elevation_sc
# Fit proposed model
width_mod3 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc + (1|genotype), data = bg_nocomp)
# Diagnostic plots
plot_model(width_mod3, type = "diag") # Looks good
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

##
## [[2]]
## [[2]]$genotype
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[3]]

##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

# See if we can add any random effects (all additive)
width_mod4 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc +
(1+co2+salinity|genotype), data = bg_nocomp) # fits
width_mod5 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc +
(1+co2+elevation_sc|genotype), data = bg_nocomp)
## boundary (singular) fit: see ?isSingular
width_mod6 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc +
(1+salinity+elevation_sc|genotype), data = bg_nocomp) # fits
width_mod7 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc +
(1+salinity|genotype), data = bg_nocomp) # fits
width_mod8 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc +
(1+elevation_sc|genotype), data = bg_nocomp) # fits
width_mod9 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc +
(1+co2|genotype), data = bg_nocomp)
## boundary (singular) fit: see ?isSingular
anova(width_mod3, width_mod4)
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod3: location * elevation_sc + (1 | genotype)
## width_mod4: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod4: location * elevation_sc + (1 + co2 + salinity | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## width_mod3 11 331.06 368.83 -154.53 309.06
## width_mod4 16 337.85 392.79 -152.93 305.86 3.2055 5 0.6683
anova(width_mod3, width_mod6)
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod3: location * elevation_sc + (1 | genotype)
## width_mod6: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod6: location * elevation_sc + (1 + salinity + elevation_sc |
## width_mod6: genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## width_mod3 11 331.06 368.83 -154.53 309.06
## width_mod6 16 336.74 391.68 -152.37 304.74 4.3164 5 0.5048
anova(width_mod3, width_mod7)
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod3: location * elevation_sc + (1 | genotype)
## width_mod7: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod7: location * elevation_sc + (1 + salinity | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## width_mod3 11 331.06 368.83 -154.53 309.06
## width_mod7 13 332.58 377.22 -153.29 306.58 2.4836 2 0.2889
anova(width_mod3, width_mod8)
## refitting model(s) with ML (instead of REML)
## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod3: location * elevation_sc + (1 | genotype)
## width_mod8: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity +
## width_mod8: location * elevation_sc + (1 + elevation_sc | genotype)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## width_mod3 11 331.06 368.83 -154.53 309.06
## width_mod8 13 333.72 378.36 -153.86 307.72 1.3375 2 0.5123
# None are significant
# Plot the fixed effects
Anova(width_mod3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: width_scam_mid
## Chisq Df Pr(>Chisq)
## date_planted_grp 8.4408 1 0.003669 **
## age 0.3273 1 0.567251
## co2 0.0144 1 0.904454
## location 0.0168 1 0.896926
## salinity 1.1383 1 0.286019
## elevation_sc 181.7487 1 < 2.2e-16 ***
## location:salinity 4.2430 1 0.039412 *
## location:elevation_sc 7.7457 1 0.005384 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Make fixed effect plot for location by elevation
width_means <- emmeans(width_mod3, ~location:elevation_sc,
at= list(elevation_sc = seq(-1.72, 1.54, 0.1)), type = "response")
width_means_out <- summary(width_means)
bg_nocomp_plot %>%
mutate(location = case_when(location == "corn" ~ "Corn Island",
T ~ "Kirkpatrick Marsh")) -> bg_nocomp_plot
width_means_out %>%
mutate(location = case_when(location == "corn" ~ "Corn Island",
T ~ "Kirkpatrick Marsh")) %>%
mutate(elevation = elevation_sc*scaling + centering) %>%
ggplot(aes(x = elevation, y = emmean, color = location)) +
geom_line(size = 1.5) +
geom_point(data = bg_nocomp_plot, aes(x = elevation, y = width_scam_mid), alpha = 0.5) +
theme_bw() + scale_color_manual(values = c("#1b9e77", "#7570b3")) +
ylab("Mean stem width (mm)") +
xlab("Elevation (m NAVD88)") +
guides(color=guide_legend("")) +
theme(legend.position = "top") -> m
# Same but for location:salinity interaction
width_mod3 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
location*salinity + location*elevation_sc + (1|genotype), data = bg_nocomp_plot)
emmip(width_mod3, location ~ salinity, CI = T, weights = "proportional") +
geom_jitter(data = bg_nocomp_rs_plot, aes(x = salinity, y = width_scam_mid, color = location),
height = 0, width = 0.1, alpha = 0.3, shape = 1) +
ylab("Mean stem width (mm)") +
scale_color_manual(values = c("#1b9e77", "#7570b3")) + theme_bw() +
theme(legend.position = "top") + guides(color=guide_legend("")) +
scale_size_manual(values=c(2,2)) -> n
# Random effects plot
REsim(width_mod3) %>%
mutate(groupID = fct_reorder(groupID, mean)) %>%
ggplot(aes(x = groupID)) +
geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
coord_flip() +
ylab("Random effect") +
xlab("Genotype") +
geom_hline(aes(yintercept = 0), linetype = "dashed") +
theme_bw() + scale_x_discrete(position = "top") +
theme(legend.position = "none") +
scale_color_manual(values = col_vector) -> o
mn <- cowplot::plot_grid(m, n, ncol = 1, labels = c("a", "b"))
## Warning: Removed 1 rows containing missing values (geom_point).
o_alone <- cowplot::plot_grid(o, labels = c("c"))
mn+o_alone+plot_layout(widths = c(4,1)) -> width_plot
# Fit a model without age, location or genotype
width_mod_null <- lm(width_scam_mid ~ date_planted_grp + co2 +
salinity + elevation_sc, data = bg_nocomp_plot)
width_null_r2 <- r.squaredGLMM(width_mod_null)[2]
width_r2 <- r.squaredGLMM(width_mod3)[2]
Belowground biomass
bgb_mod <- lmer(total_bg ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|site_frame) + (1|genotype), data = bg_nocomp)
#step(bgb_mod)
VarCorr(bgb_mod)
## Groups Name Std.Dev.
## genotype (Intercept) 1.40014
## site_frame (Intercept) 0.27376
## Residual 3.38998
bgb_mod2 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + origin_lab +
(age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
(1|genotype), data = bg_nocomp)
step(bgb_mod2, keep = c("co2", "elevation_sc", "salinity", "age", "location"))
## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 38 -565.32 1206.7
## (1 | genotype) 0 37 -568.81 1211.6 6.9719 1 0.00828 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF
## age:location:co2:salinity:elevation_sc 1 0.30 0.30 1 191.632
## origin_lab 2 0.80 0.80 1 31.107
## age:location:co2:elevation_sc 3 2.13 2.13 1 183.313
## age:location:co2:salinity 4 6.42 6.42 1 176.485
## age:location:salinity:elevation_sc 5 8.66 8.66 1 180.800
## age:location:elevation_sc 6 0.00 0.00 1 181.213
## age:location:salinity 7 3.31 3.31 1 185.205
## age:co2:salinity:elevation_sc 8 13.57 13.57 1 196.165
## age:co2:elevation_sc 9 13.05 13.05 1 190.620
## location:co2:salinity:elevation_sc 10 13.26 13.26 1 200.299
## co2:salinity:elevation_sc 11 0.65 0.65 1 200.325
## location:co2:elevation_sc 12 2.49 2.49 1 194.922
## location:co2:salinity 13 2.76 2.76 1 185.776
## age:salinity:elevation_sc 14 11.66 11.66 1 190.162
## age:elevation_sc 15 6.95 6.95 1 190.408
## location:salinity:elevation_sc 16 11.15 11.15 1 193.654
## location:elevation_sc 17 0.10 0.10 1 191.214
## location:salinity 18 2.30 2.30 1 196.923
## co2:elevation_sc 19 12.25 12.25 1 200.383
## salinity:elevation_sc 20 24.80 24.80 1 197.169
## age:location:co2 21 24.55 24.55 1 202.999
## location:co2 22 1.64 1.64 1 204.121
## age:location 23 11.33 11.33 1 25.241
## weight_init_sc 0 53.05 53.05 1 213.800
## date_planted_grp 0 53.89 53.89 1 205.811
## location 0 1.50 1.50 1 27.270
## elevation_sc 0 0.02 0.02 1 199.343
## I(elevation_sc^2) 0 386.08 386.08 1 204.957
## age:co2:salinity 0 45.46 45.46 1 195.963
## F value Pr(>F)
## age:location:co2:salinity:elevation_sc 0.0260 0.87200
## origin_lab 0.0699 0.79323
## age:location:co2:elevation_sc 0.1854 0.66726
## age:location:co2:salinity 0.5606 0.45500
## age:location:salinity:elevation_sc 0.7585 0.38496
## age:location:elevation_sc 0.0000 0.99442
## age:location:salinity 0.2915 0.58994
## age:co2:salinity:elevation_sc 1.1995 0.27476
## age:co2:elevation_sc 1.1492 0.28508
## location:co2:salinity:elevation_sc 1.1603 0.28270
## co2:salinity:elevation_sc 0.0563 0.81261
## location:co2:elevation_sc 0.2169 0.64193
## location:co2:salinity 0.2425 0.62296
## age:salinity:elevation_sc 1.0274 0.31206
## age:elevation_sc 0.6120 0.43500
## location:salinity:elevation_sc 0.9874 0.32161
## location:elevation_sc 0.0093 0.92346
## location:salinity 0.2048 0.65138
## co2:elevation_sc 1.0961 0.29637
## salinity:elevation_sc 2.2101 0.13871
## age:location:co2 2.1710 0.14218
## location:co2 0.1442 0.70451
## age:location 0.9984 0.32718
## weight_init_sc 4.6787 0.03165 *
## date_planted_grp 4.7525 0.03039 *
## location 0.1321 0.71908
## elevation_sc 0.0021 0.96312
## I(elevation_sc^2) 34.0485 2.075e-08 ***
## age:co2:salinity 4.0094 0.04663 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## total_bg ~ weight_init_sc + date_planted_grp + age + location +
## co2 + salinity + elevation_sc + I(elevation_sc^2) + (1 |
## genotype) + age:co2 + age:salinity + co2:salinity + age:co2:salinity
bgb_mod3 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation +
+ location + I(elevation^2) + (1 | genotype) + age*co2*salinity, data = bg_nocomp_plot)
# Diagnostics
plot_model(bgb_mod3, type = "diag") # Pretty!
## [[1]]
## `geom_smooth()` using formula 'y ~ x'

##
## [[2]]
## [[2]]$genotype
## `geom_smooth()` using formula 'y ~ x'

##
##
## [[3]]

##
## [[4]]
## `geom_smooth()` using formula 'y ~ x'

# Fixed effects plots
car::Anova(bgb_mod3)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: total_bg
## Chisq Df Pr(>Chisq)
## weight_init_sc 4.6787 1 0.030539 *
## date_planted_grp 4.7525 1 0.029256 *
## age 0.0091 1 0.923897
## co2 13.8404 1 0.000199 ***
## salinity 2.6561 1 0.103153
## elevation 32.5434 1 1.166e-08 ***
## location 0.1321 1 0.716282
## I(elevation^2) 34.0485 1 5.376e-09 ***
## age:co2 0.0088 1 0.925363
## age:salinity 0.2140 1 0.643665
## co2:salinity 2.8146 1 0.093409 .
## age:co2:salinity 4.0094 1 0.045247 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bgb_means <- emmeans::emmeans(bgb_mod3, ~co2:salinity:age:elevation, at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")
bgb_means_out <- summary(bgb_means)
bgb_means_out %>%
ggplot(aes(x = elevation, y = emmean, color = co2)) +
geom_line(size = 1.5, alpha= 0.8) +
geom_line(aes(x = elevation, y = lower.CL), linetype = "dashed") +
geom_line(aes(x = elevation, y = upper.CL), linetype = "dashed") +
facet_grid(age~salinity) +
geom_point(data = bg_nocomp_plot, aes(x = elevation, y = total_bg), alpha = 0.5) +
theme_bw() +
scale_color_manual(values = c("gray47", "#66a61e")) +
ylab("Belowground biomass (g)") +
xlab("Elevation (m NAVD88)") +
guides(color=guide_legend(title=expression(CO[2]))) +
theme(legend.position = "top") -> p
# See if we can add random effects
bgb_mod4 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + salinity + elevation_sc | genotype) + age*co2*salinity, data = bg_nocomp_plot)
## boundary (singular) fit: see ?isSingular
bgb_mod5 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + co2 + elevation_sc | genotype) + age*co2*salinity, data = bg_nocomp_plot)
## boundary (singular) fit: see ?isSingular
bgb_mod6 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + salinity + co2 | genotype) + age*co2*salinity, data = bg_nocomp_plot)
## boundary (singular) fit: see ?isSingular
bgb_mod7 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + co2 | genotype) + age*co2*salinity, data = bg_nocomp_plot)
## boundary (singular) fit: see ?isSingular
bgb_mod8 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + salinity | genotype) + age*co2*salinity, data = bg_nocomp_plot)
## boundary (singular) fit: see ?isSingular
bgb_mod9 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + elevation_sc | genotype) + age*co2*salinity, data = bg_nocomp_plot)
## boundary (singular) fit: see ?isSingular
# No random slopes
# Create caterpillar graph
REsim(bgb_mod3) %>%
mutate(groupID = fct_reorder(groupID, mean)) %>%
ggplot(aes(x = groupID)) +
geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
coord_flip() +
ylab("Random effect") +
xlab("Genotype") +
geom_hline(aes(yintercept = 0), linetype = "dashed") +
theme_bw() + scale_x_discrete(position = "top") +
theme(legend.position = "none") +
scale_color_manual(values = col_vector)-> q
p + q + plot_layout(widths = c(4,1)) -> bgb_plot
# Build null model
bgb_null <- lm(total_bg ~ weight_init_sc + date_planted_grp + elevation +
+ I(elevation^2) + co2*salinity, data = bg_nocomp_plot)
bg_null_r2 <- r.squaredGLMM(bgb_null)[2]
bg_r2 <- r.squaredGLMM(bgb_mod3)[2]


## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

|
trait
|
null r2
|
evo r2
|
|
aboveground biomass
|
0.35
|
0.46
|
|
root-to-shoot ratio
|
0.60
|
0.69
|
|
belowground biomass
|
0.23
|
0.32
|
|
root distribution parameter
|
0.46
|
0.76
|
|
mean stem height
|
0.38
|
0.60
|
|
mean stem width
|
0.42
|
0.62
|